@jngrad what is the difference between the tests https://github.com/espressomd/espresso/blob/python/testsuite/python/constant_pH.py and https://github.com/espressomd/espresso/blob/python/testsuite/python/constant_pH_stats.py ? Can one be removed?

The former checks the interface using a chemical reaction with cherry-picked parameters to run in a few seconds in coverage builds. The latter is a statistical test that checks convergence. They are otherwise almost identical, and constant_pH_stats.py will probably be removed in the future.

@jngrad I get a message about some cpH core tests failing which you introduced, could you please have a look? I guess the values are hardcored and do not match anymore.

@jngrad I get a message about some cpH core tests failing which you introduced, could you please have a look? I guess the values are hardcored and do not match anymore.

You will need to adapt the acceptance rate in the unit test with the result of auto const g = calculate_factorial_expression_cpH(...).

But before you spend more time on this, the core team should first decide whether we really want to introduce this silent breaking change to the Constant pH algorithm.

The statistical properties of the algorithm before and after the change are shown to be the same in the referred paper. I would be in favour of delivering this change.

@jonaslandsgesell can you please explain how the asymmetric probabilities fix the bug when the reaction is not possible in either direction, as described in #4197 As far as I understand, moving from symmetric to asymmetric proposal probabilities still ensures that the sum of probabilities is equal to 1.0.

The statistical properties of the algorithm before and after the change are shown to be the same in the referred paper. I would be in favour of delivering this change.

I am convinced that it works on paper but I think that the implementation should be tested on various non-trivial cases before it is merged. I am afraid that the standard test cases might not be sufficient for that. We should also verify its impact on performance.

@RudolfWeeber proposed to have a discussion on the future of reaction ensemble. This is a big change that affects the core of the reaction ensemble, so I suggest to wait with this merge until this discussion takes place.

@kosovan all_reactant_particles_exist() is checked in generic_oneway_reaction. If it returns false no reaction is attempted because a reaction step removing the non-existing reactant particles is not possible. You can also see this via constructing any situation you like and attaching a debugger. The problem in the previous code was the way the asymmetric proposal probability was implemented.

I would not defer this merge because it
A) removes some bogous code which made the proposal probabilities asymmetric.
B) If you want to make other changes please do not introduce a coupling with regards to this PR. The changes you want to make should be regarded independently.

What non-trivial test you have in mind?

A) removes some bogous code which made the proposal probabilities asymmetric.

I am not objecting against these changes, I just want to be sure that they are thoroughly tested.

B) If you want to make other changes please do not introduce a coupling with regards to this PR. The changes you want to make should be regarded independently.

I do not understand. If I just mention another issue, does github automatically treat them as dependent?

What non-trivial test you have in mind?

I need some time to think of that. I was not aware that these changes were planned, so I have not thought of it before.

Some points that come to my mind immediately: The standard test cases are not exactly reproducible, they have to rely on some arbitrarily chosen tolerance. This change affects the acceptance probability, and if there is a bug, it might not be immediately visible with a tolerance that is set so that the test completes quickly. Therefore, I would try to generate high-quality data using the old and the new implementation, and check that they produce identical results. Definitely, I would focus on various corner cases, with alpha approx. 0 and alpha approx. 1, and I would like to test if/how the performance has changed with the new code.

Some points that come to my mind immediately: The standard test cases are not exactly reproducible, they have to rely on some arbitrarily chosen tolerance. This change affects the acceptance probability, and if there is a bug, it might not be immediately visible with a tolerance that is set so that the test completes quickly. Therefore, I would try to generate high-quality data using the old and the new implementation, and check that they produce identical results. Definitely, I would focus on various corner cases, with alpha approx. 0 and alpha approx. 1, and I would like to test if/how the performance has changed with the new code.

If you don't trust the current test setup it's impossible to safely change any part of the reaction ensemble...

If you don't trust the current test setup it's impossible to safely change any part of the reaction ensemble...

The current cpH test completes on my laptop in 2.807s after it runs 1000 samples for one specific set of conditions. As far as I remember, it has been chosen as a compromise between accuracy and speed, with a preference for speed. This is sufficient in many cases, if you want to check that the code does not crash due to a programming error but it is not sufficient to check that the physics has not changed.

I think that one needs different levels of cautiousness, depending on the possible impact a particular change has. The change in this pull request has similar impact on the cpH calculations as if you would replace the integrator scheme or a random number generator. It is definitely worth more testing before it comes to production.

For example, the test case only tests an ideal system (no interactions), so it cannot catch even such a simple mistake as a wrong sign before Delta U in the acceptance probability.

For example, the test case only tests an ideal system (no interactions), so it cannot catch even such a simple mistake as a wrong sign before Delta U in the acceptance probability.

I know you just made an example, but the changes introduced do not affect interactions because the factor involving interactions is untouched...

I know you just made an example, but the changes introduced do not affect interactions because the factor involving interactions is untouched...

It should not affect interactions. But maybe it affects something else that you have not considered. Just give us some time to test it on a non-trivial problem.

Just to be sure, I can perform simulations for some interacting systems in some non-trivial conditions and check that the results I get with your new implementation are the same that the ones I was getting with the actual code. Also could be a good chance to check the performance :)

I did some calculations of the titration curve of a weak polyacid (with WCA and electrostatic interactions). Attached you can find the results with the new implementation (label "new") and the current development code (label "devel"). In all the pH range both algorithms yield to the same average degree of ionization (\alpha) within the estimated error. (Note: the results are rounded to the third decimal).
I also tested that this new implementation does indeed solve the infinite loop bug that I mentioned in #4197, by returning a no-op. I also noticed that the new implementation seems to be slightly faster than the current develop version.
In my opinion, the branch can be safely merged.

PS: The cpp cpH unittest also checks the impact of the input energy.